library(plm)
library(lmtest)

countries = c('indo', 'colo')
dvs = c('any', 'high', 'spike')
shocks <- list(indo = c(A='rhs_tmp_RainDevL0',
                        B='rhs_tmp_TempDevL0',
                        C='rhs_shk_food_wtDln_priceL0',
                        D='rhs_shk_majcash_wtDln_priceL0'),
               colo = c(A='rhs_rain_deviationL0',
                        B='rhs_temp_deviationL0',
                        C='rhs_cofintxlinternalpL0',
                        D='rhs_oilprod88xlopL0'))

shock_check <- array(NA,
                     dim = c(length(countries),
                             length(dvs),
                             4,
                             2,
                             3),
                     dimnames = list(countries,
                                     dvs,
                                     c('A', 'B', 'C', 'D'),
                                     c('current', 'future'),
                                     c('b', 'se', 'ast')))

for (country in countries) {
  for (v in dvs) {
    for (s in c('A', 'B', 'C', 'D')) {
      print(paste0(country, ' -- ', v , ' -- ', s))
      source(paste("config_code/",country,"_data_setup.R",sep = ""))
      source(paste("config_code/",country,"_dependent_vars.R",sep = ""))
        
      X <- shocks[[country]][s]
      
      dta$std_shock = (dta[,X] - mean(dta[,X]))/sd(dta[,X])
      mod <- plm(paste0(dv,
                      '~',
                      'std_shock'),
               data=pdata.frame(dta,
                                index = c(loc.var, t.var), 
                                drop.index = F, row.names = F),
               model="within", effect="twoways")
    
      G <- length(unique(dta[,loc.var]))
      N <- length(dta[,loc.var])
      dfa <- (G/(G - 1)) * (N - 1)/mod$df.residual
      firm_c_vcov <- dfa * vcovHC(mod, type = "HC0", cluster = "group", adjust = T)
      se1 <- coeftest(mod, vcov = firm_c_vcov)
      shock_check[country, v, s,
                  'future', ] <- se1['std_shock',c('Estimate',
                                                   'Std. Error',
                                                   'Pr(>|t|)')]
      
      remove(mod)
      remove(firm_c_vcov)
      remove(se1)
      
      mod <- plm(paste0(lagdv,
                       '~',
                       'std_shock'),
                 data=pdata.frame(dta,
                             index = c(loc.var,
                                       t.var), 
                             drop.index = F, row.names = F),
                 model="within", effect="twoways")
    
      
      G <- length(unique(dta[,loc.var]))
      N <- length(dta[,loc.var])
      dfa <- (G/(G - 1)) * (N - 1)/mod$df.residual
      firm_c_vcov <- dfa * vcovHC(mod, type = "HC0", cluster = "group", adjust = T)
      se1 <- coeftest(mod, vcov = firm_c_vcov)
      shock_check[country, v, s,
                  'current', ] <- se1['std_shock',c('Estimate',
                                                    'Std. Error',
                                                    'Pr(>|t|)')]
      
      remove(mod)
      remove(firm_c_vcov)
      remove(se1)
    }
  }
}
shocks_labels <- list(indo = c(A='Rainfall',
                        B='Temperature',
                        C='Food Price',
                        D='Major Cash Crop Price'),
               colo = c(A='Rainfall',
                        B='Temperature',
                        C='Coffee Price',
                        D='Oil Price'))

outcome.labels <- c(any="(a) Any violent event",
                    high="(b) $\\ge$ 5 violent events",
                    spike="(c) $\\ge$ 1 s.d. increase in events")

dimnames(shock_check)[[2]] <- outcome.labels[dimnames(shock_check)[[2]]]

printC44(table = shock_check['indo',,,,],
         filepath = "tables",
         file = "table_C44",
         shock_key=shocks_labels[['indo']])

printC44(table = shock_check['colo',,,,],
         filepath = "tables",
         file = "table_C45",
         shock_key=shocks_labels[['colo']])

